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ABSTRACT 

A gravitational lens model is presented for the newly discovered 10-image system 
B1933+503. The underlying object, revealed by modeling, is a triple radio source on 
the scale of a couple of hundred mas that is well-aligned along the line of sight with 
a foreground and somewhat flattened lensing galaxy, whose orientation and location 
match that of an observed galaxy, known to be at a redshift of 0.755. Uncertain- 
ties in the modeling are obtained by a Monte Carlo exercise. Observational tests of 
the lens model are proposed, and the time delays between various pairs of images 
are determined as the core of the source is known to be significantly variable. Fu- 
ture observations of the lens hold the key to using B1933+503 to constrain Hubble's 
Constant. Despite the absence of a source redshift, the system's utility as a probe of 
the lens galaxy's structure is unparalleled as it provides a surfeit of easily identifiable 
constraints for modeling the system. 
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1 INTRODUCTION 

In a companion paper, Sykes et al. (1997) report the discov- 
ery of a spectacular arcsecond-scale gravitationally lensed 
radio system, B1933+503 (1934+504) in J2000 coordinates), 
which was found in the course of the Cosmic Lens All-Sky 
Survey (CLASS). CLASS is a search for radio sources that 
exhibit multiple components with flat radio spectral indices 
(a < 0.5, with flux density S v ~ v~ a ). Such systems are typ- 
ically extragalactic, and result from gravitational multiple 
imaging of a background source such as a quasar or an AGN 
by a foreground galaxy lying along the line of sight to it. The 
first gravitational lens to be identified on the basis of its mul- 
tiple flat-spectrum radio components was PKS1830— 211, by 
Rao & Subrahmanyan (1988). A survey conducted from Jo- 
drell Bank (Patnaik et al. 1992), now known as the Jodrell 
VLA Astrometric Survey or JVAS, employed this property 
as a highly successful filter for identifying new lensed sys- 
tems: B0218+357 (Patnaik et al. 1993), B1422+231 (Pat- 
naik et al. 1992), B1030+074 (King & Browne 1996 and 
Xanthopoulos et al. 1997) and B1938+666 (King et al. 
1997a, King et al. 1997b), in addition to the promising can- 
didate lensed system B2114+022 (Augusto et al. 1996). The 
ongoing CLASS project follows on from and contains JVAS, 
and together these have discovered more than 13 systems 
so far (inclusive of candidates yet to be confirmed) in a 
survey of nearly 8000 sources. Descriptions of this survey 
are to be found in Myers et al.(1995), Jackson et al. (1995, 
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1997b) and Browne et al. (1997). Newly established finds in- 
clude B1608+656 (Myers et al. 1995, Fassnacht et al. 1996), 
B1600+434 (Jackson et al. 1995, Koopmans et al. 1997a), 
B0712+472 (Jackson et al. 1997a), B1933+503 (Sykes et al. 
1997) and B1130+382 (Koopmans et al. 1997b). 

B1933+503 bears a superficial resemblance to an Ein- 
stein Ring in its morphology, but is composed of at least 
10 discrete components, as reported by Sykes et al. (1997). 
It has been investigated over the range of radio frequencies 
1.7-15 GHz (VLA 5, 8.4 & 15 GHz, MERLIN 1.7 & 5 GHz 
and VLBA 1.7 & 5 GHz); only the MERLIN 1.7 GHz obser- 
vations appear to have the combination of dynamic range 
and resolution required to pick up all 10 components. How- 
ever, Sykes et al. (1997) provide sufficient observational ev- 
idence to note that (a) four of the components have similar, 
complex radio spectra peaking around 5 GHz and are com- 
pact, being a few mas in size, (b) of the remaining six, four 
have steep radio spectra over the range of frequencies ob- 
served, but appear to be fairly compact (also on the scale of 
a few mas). The remaining two components probably also 
have steep radio spectra. Sykes et al. (1997) also infer that 
the components which have complex radio spectra appear to 
be variable by as much as 33% at 15 GHz over the timescale 
of a couple of months. HST observations with WFPC2 suc- 
ceed in picking up what appears to be a flattened galaxy 
located near the centre of the ring of radio components. 
This object has an I magnitude of 20.6 at 810 nm and is 
elongated at a position angle of —40 ± 5 degrees, with an 
axial ratio (b/a) in the plane of the sky of 0.45 — 0.55. At 
540 nm, there is no sign of the galaxy down to V magnitude 
22.5. An optical spectrum taken with the Keck telescope us- 
ing the Low Resolution Imaging Spectrometer (LRIS) shows 
the presence of absorption and emission corresponding to a 
redshift of 0.755, which probably corresponds to the galaxy 
rather than the source since no optical counterparts of the 
images have been detected as yet (Sykes et al. 1997). 



2 B1933+503: THE LENS INTERPRETATION 

A galaxy acting as a lens typically produces five, three or 
a single image(s) of a background source, with increasing 
degree of misalignment of the source in the plane of the sky 
from the line of sight to the lens (see, e.g., Schneider, Ehlers 
& Falco 1992 for a review of the properties of gravitational 
lensing by galaxies). When five or three images are formed, 
one occurs very near the centre of the galaxy and suffers 
a high degree of demagnification, so observationally there 
would appear to be either four images (a quad) or two im- 
ages (a double) of the source. Most radio sources discovered 
so far that are lensed by galaxies and exhibit multiple com- 
pact flat-spectrum features can be clearly classified as either 
quad or double systems. B1938+666 (King et al. 1997) has 
two components in the background source, one of which is 
quadruply imaged while the other is a double. With its 10 
features, B1933+503 is thus simplest understood as a triple 
radio source, the individual components of which have been 
multiply imaged into a quad, a quad and a double. 

The MERLIN 1.7 GHz map of Sykes ct al. (1997) (at 
top left in fig.l of that paper) is the starting point of the 
present modeling exercise. Most of the components in this 
map can be grouped into one or the other of two quad con- 



figurations. The first quad involves features 2, 5 and 7. The 
elongated morphology of feature 2 may be interpreted as a 
pair of images that are partially merged across the tangen- 
tial critical curve (see Fig. 1(a)). There appear to be two 
strong peaks of flux density in feature 2 in the 1.7 GHz map 
of Sykes et al. (1997), and these will be referred to as 2a 
(east) and 2b (west) in the present work. The second quad 
is formed by components 1, 3, 4 and 6. Component 8, lying 
as it does within the circle of images, cannot then be singly 
imaged; its counterpart image is to be found in the faint 
feature la to the NE of component 1. 

The two-image configuration formed by la and 8 is un- 
usual in that 8 is so bright by comparison to la, despite its 
being nearer to the lens centre. This provides an interesting 
constraint on the lens model, since this relatively rare con- 
figuration can be obtained if image 8 is derived from three 
images that have just barely merged with each other into 
a single bright image. The source that is imaged into la 
and 8 must lie just outside a cusp of the tangential caus- 
tic in the source plane, but within the radial caustic. (This 
is illustrated in Fig. 1(b), a result from the lens modeling 
described in Sects. 3 and 4). 

Within this picture, the underlying source prior to 
imaging by the lens would consist of a central radio core 
which shares the spectral and morphological properties of 
images 1, 3, 4, and 6, and is thus flat-spectrum and compact. 
The radio core is flanked by two steeper spectrum radio fea- 
tures, the trio being found by actual modeling to lie almost 
in a straight line, on the scale of a couple of hundred mas 
(Fig. 1(b)). 

In principle, an alternative classification of the images 
could group 2 (seen again as a pair of images), 7, and 4 
together, and 1, 3, 6, and 8 into a second quad. This would 
require 5 to be singly-imaged, as also la, unless each of them 
is an element of a double with a counterpart image within 
the circle of images which is too weak to be mapped. We 
can discard this and related scenarios in view of the support 
that the earlier classification receives from the radio spectra 
in figure 2 of Sykes et al. 1997, the VLBA 5GHz map in 
figure 1 (Sykes et al. 1997), and the pattern of variability 
discussed in that paper. (Amusingly, the image classification 
described earlier was worked out prior to any knowledge of 
the radio spectra or the VLBA observations.) 



3 LENS MODELING 

A modeling code has been developed that seeks a best- 
fit model for B1933+503, following the methods of Kayser 
& Schramm (1988) and Kochanek (1991). These employ a 
penalty function, which is minimised over the parameter 
space of a parametrized lens mass model, to yield a 'best 
fit' lens. Although the observed flattening of the lens sug- 
gests that it could be a disk galaxy which would typically 
consist of a disk, a spheroid and a halo, it is treated here as 
a single component elliptical lens approximating the overall 
mass distribution. The lensing galaxy is described by a non- 
singular isothermal ellipsoidal mass profile (the 'PIEMD' of 
Kassiola & Kovner 1993). The mass density distribution in 
the lens, p, follows the form: 

p(m) = Po /{\ + (m/af), (1) 
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where a is the scale length of the mass distribution, p a is 
the central mass density, and m 2 — x 2 + y 2 / (I — e 2 ) + z 2 is 
the semi-major axis of an ellipsoidal shell of constant mass 
density. The lens is thus an oblate spheroid with axial ratio 
given by \/I — e 2 , and is assumed to be viewed edge-on. (The 
x— and y— axes are confined to the plane of the sky, the 
orientation being set by modeling the lensed system; in the 
absence of concrete information about the nature of the lens, 
this model approximates the behaviour of either an elliptical 
or a spiral galaxy). The lens model has six parameters, these 
being the coordinates of the lens centre in the lens plane, a 
mass parameter describing the strength of the lens, its scale 
length, ellipticity and orientation in the plane of the sky. 
The mass parameter, a m , is related to the central density 
and the scale length: p = 9<r 2 n /47rGa 2 . 

The complex scattering function formalism of Bourassa, 
Kantowski & Norton (1973) and Bourassa & Kantowski 
(1975) is employed, which permits analytical expressions to 
be obtained for the lensing action of a spheroidal mass dis- 
tribution. In this formalism, the lens equation relating the 
source position in the plane of the sky, z s = x s + iy s , to that 
of its images, z; = x + iy, via a scattering function, I(x,y), 
is given by: 



I* denotes the complex conjugate of the scattering func- 
tion, which is proportional to aj,. I* is a function of image 
position Zi as well as of the lens parameters pj. G is the 
gravitational constant, and Fd is the ratio of two angular 
diameter distances, that between lens and source, Di s , and 
that between observer and source, D s (since the redshift of 
the source is unknown at present, F d is absorbed into the 
lens parameter cr^ and the image and source positions are in 
angular units in the plane of the sky). The scattering func- 
tion for the density profile in Eqn.(l) is given by expression 
(4.1.2) of Kassiola & Kovner (1993). 

3.1 Constraints on the Modeling 

Details of the image positions, relative to an initial guess 
lens centre, are listed in Table 1, in which the images are 
grouped according to their membership of a particular im- 
age configuration. Each multiple image configuration of N 
images provides 2(N-1) positional constraints (eliminating 
the common source position between pairs of images in a 
group; see Eqn.(2)). Thus, the two quads and the double 
image configuration of B1933+503 supply a maximum of 
6+6+2=14 constraints. This surfeit of constraints over pa- 
rameters permits the luxury of ignoring the relatively un- 
certain positions of images 2a and 2b: these images straddle 
the tangential critical curve and each should theoretically 
be significantly brighter than either 5 or 7, which share the 
same source but are removed from any critical curves. This 
source, being of steep radio spectrum, is presumably non- 
variable (note the constancy of the flux ratio for images 7 
and 5 in Table 1; the observations have been taken at differ- 
ent epochs). However, from table 2 of Sykes et al. (1997) it 
is apparent that in almost every observation, the combined 
flux density of 2a and 2b is lower than that of 5 or 7. One is 
led to conclude that the source of the image set (2a, 2b, 5, 
7) is extended and is only partially imaged in the vicinity of 



the tangential critical curve. Thus both the positions of 2a 
and 2b and their flux densities are undependable inputs to 
the modeling. Nevertheless, the fact that the image parities 
must be reversed between the positions of features 2a and 2b 
(see Table 2) proves very useful in constraining the location 
of the tangential critical curve during the modeling process, 
and this parity information is usable. Thus, only 2+6+2=10 
positional constraints are actually employed. 

Ratios between the observed flux densities of pairs of 
images in a particular configuration generally provide an 
effective set of constraints on the modeling. However, if 
the source shows temporal variations in flux density, such 
changes manifest at different times in the various images 
owing to differential time delays in the arrival time at the 
observer for light from them. In B1933+503, it is known 
that the quad images (1, 3, 4, 6), arising from the flat- 
spectrum core, are significantly variable as discussed in Sec- 
tion 1. Hence constraints derived from their flux densities at 
a given epoch of observation could be uncertain by as much 
as 40% in the higher frequency observations. Accordingly, 
these are not be used in the present modeling exercise. Of 
the quad (2a, 2b, 5, 7), flux density ratios involving 2a and 
2b are neglected, as described earlier, but the ratio for im- 
age 7 to image 5 appears to be robust (see Table 1); the 
source is of steep-spectrum and is not expected to be vari- 
able. Similarly, the flux density ratio of image 8 to la should 
be non-variable, but in practice this is affected particularly 
by the uncertainty in la, on account of its faintness. No 
flux density ratio constraints are actually employed in the 
modeling, but a good model may reasonably be expected to 
reproduce the ratios 7/5 and 8/la. 

In the absence of reliable image flux density ratios, the 
image parities provide important constraints on the model- 
ing. The parity ( — 1)™ for an image in a given configuration 
(quad or double) represents how the source maps onto that 
image, n being the number of reflections that it experiences 
about a set of coordinate axes centred on the image (for a 
single lens plane, no rotations occur). The theoretically ex- 
pected image parities can be easily worked out by inspection 
for standard elliptical lenses (see Blandford & Narayan 1986, 
for example), and are listed in Table 2 for B1933+503. This 
set of constraints essentially determines the lens orientation 
in the plane of the sky. 

3.2 The Penalty Function 

For the desired lens model, the positions of those images be- 
longing to the same quad or double configuration must map 
back through the lens, via the lens equation, to their shared 
source position. The penalty function to be minimised is ob- 
tained in the following manner : In general, having guessed 
an initial set of lens parameters, image positions correspond- 
ing to the same quad or double do not map back via Eqn.(2) 
to the same source position. A summed squared mismatch 
between the various source positions obtained by backpro- 
jecting the observed image positions, taken pairwise and di- 
vided by the number of images involved, forms the first part 
of the penalty function; this should ideally be as small as 
possible for each set of quads and the double. In the liter- 
ature, it has been popular to multiply the errors in source 
recovery (departures from the average recovered source posi- 
tion) by the model magnifications of each image, converting 
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these errors to values that can be compared directly with the 
observational errors (e.g. Kochanek 1991). The minimization 
of the error function is then carried out in the image plane. 
In the present work, we have avoided using the model mag- 
nifications to bias the penalty function. This is with good 
reason: trials in which it was attempted to minimize the mis- 
match in the image plane after multiplying the source plane 
mismatches with the local image magnifications tended to 
locate models with enormously high image magnifications 
(~ 10 4 and greater). In fact, given the freedom that a mass 
model with a finite scale length permits, the program tended 
to place the tangential critical curve as close to the images 
as it could (by centering the lens with respect to the images 
and making it as round as possible) . Trying to minimize the 
source plane mismatch appears to be free of this bias and the 
results accord with trials with synthetic data from (known) 
lens models, obtained using a root-finding algorithm to solve 
the lens equation for the image positions corresponding to 
a given source. 

The second part of the penalty function incorporates 
the image parity constraints. Only those lens models are 
accepted which exactly match the predicted pattern of image 
parities. This ensures that unphysical configurations are not 
selected by the program. Accordingly, the penalty associated 
with a parity mismatch is very high. The penalty function 
reads as: 



E 



i Nk 

iVfc E 1 

•<J 
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+ io 4 r(P m p ei 

m— 1 



(3) 



where Nk denotes the number of images in the image group 
(quad or double) , a is the tolerance in recovering the source 
position, and P m and P em are the model and expected val- 
ues of parity for image n in a given group. The parities take 
values of either +1 or — 1. F(P m P em ) is a function with value 
zero for P m P em > 0, else it is unity. The large penalty as- 
sociated with this part of the function <f> ensures that parity 
match is obtained almost at the outset of the minimization 
process. 

Minimization of the penalty function is achieved 



using a program, 'simann.f, written by 
on Goffe et al. 1994, and available from 
public access site on the Inter net at 
http:/ /www.netlib.org/opt/simann.l). This 



Goffe (based 
the NETLIB 
the location 
program em- 
ploys the method of simulated annealing for optimization. It 
is found to be particularly useful in accomodating the sud- 
den discontinuities that arise in the penalty function, accom- 
panying changes in the model values of image parity as the 
six-dimensional parameter space is searched. Although com- 
putationally more expensive than the more popular simplex 
method § based optimization algorithms, it is significantly 
less prone to getting trapped in local optima. 



S Descriptions of both simulated annealing and simplex methods 
in optimization routines are to be found in Press et al. (1992) 



4 RESULTS AND A DISCUSSION OF THE 
LENS MODEL 

The source positions, as found for the best-fit (cf> = 9.83 
with a target a of 5 mas) lens model are displayed in Table 
2, and are plotted in Figure 1(b). Note particularly the suc- 
cess of the model in recovering the lobe double's source to 
lie just outside the cusp of the tangential critical curve, as 
mentioned in Section 2. The model image magnifications are 
also given in Table 2. The parameters of the best-fit model 
are given below; the errors quoted in brackets alongside are 
90% confidence intervals from a Monte Carlo exercise. For 
this, the image positions have been subjected to random 
deviations, the magnitudes of which follow a normal proba- 
bility distribution with a standard deviation of 5 mas from 
their nominal values. In the present case, 10000 such image 
configurations were sampled and put through the modeling 
process described above, with all six lens parameters free 
to vary. This took a total of 65 hours cpu time on a Sparc 
Ultra 1 machine. While 10 4 trials may seem a sparse sam- 
pling of the possible configuration space of random devia- 
tions in all the image positions put together, in practice the 
estimated confidence intervals alter only marginally between 
experiments with 10 4 trials and half that number (though a 
significant change occurs between 500 and 5000 trials). 

The best-fit model obtained by the methods described 
in the previous sections yields lens parameters as below 
(with 90% confidence intervals quoted alongside): 

Scale Length a: 113 mas (59, 224)mas 

Eccentricity e: 0.81 (0.73, 0.84) 

Mass Parameter a m : 79.8 (78.2,84.1) km/s 

P.A. of Lens Major Axis: -46.5° (-47.0°, -46.1°) 

Lens Centre (RA, Dec. wrt cpt. 4): (423±| , 270lg) mas 



Histograms of the distributions of the mismatch func- 
tion (j> an d the various lens parameters, as obtained from 
the Monte Carlo exercise, are presented in Figure 2. The 
'shoulder' in the plot for the mismatch function seen around 
the values of 12 to 13 is probably related to a similar fea- 
ture seen in the plot for the position angle around —46.4 
to —46.3 degrees, which is likely to be a consequence of the 
image parity constraints (sudden discontinuities result when 
a change of parameters places an image in a region of image 
space with the wrong parity). 

The circular velocity V C (R) in the equatorial plane of 
the oblate spheroidal lens galaxy with mass density distri- 
bution as in Eqn.(l), can be estimated as a function of radial 
coordinate R from expression (2-91) in Binney & Tremaine 
(1987). 



V C 2 (R) = 9o-; 



VT- 



tan 



-ivT 



V 1 + W a f 



(4) 



v /l + P 2 /(ae) 2 

In the limit of e — > 0, the right hand side of 
Eqn.(4) reduces to 9a 2 ra (D s /A s ){l - (a/P)tan _1 (P/a)}. 
Asymptotically, the right hand side of Eqn.(4) is: 
9ol(VT^/e)(D s /D ls ){n/2 - t&n- 1 {^1^ / e)} . From 
the above, it is seen that the circular velocity at a radius of 
1" from the lens centre is 183 yj D s /Di s km/s, which, as R 
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tends to large values, achieves a maximum of 198^/D S /A S 
km/s. 

The eccentricity e translates to a flattening or b/a ratio 
in the lens plane of 0.59lg'Qg. This is marginally rounder 
than the value of 0.45 — 0.55 suggested by the light pro- 
file in the HST I-band image, which is not surprising if 
the degree of flattening of the non-luminous component of 
the lens galaxy is somewhat smaller than that of its disk. 
It is of interest that the orientation of the lensing galaxy 
(p. a. —40 ± 5 degrees) accords so well with that required by 
our lens model. This strongly suggests that B1933+503 is 
lensed in the main by an isolated disk system. From Table 2, 
a comparison of the image flux density ratios with the model 
magnification values shows that the model actually predicts 
quite reasonably the flux density ratio of components 7 and 
5, and with a lesser degree of success, the somewhat less 
reliable flux density ratio 8/la. In the case of the core im- 
ages, departures of the model magnification ratios from the 
observed flux density ratios are apparent; these are proba- 
bly the result of rapid intrinsic variations in the core flux 
density, as discussed earlier. 



5 PREDICTED TIME DELAYS FOR THE 
IMAGES 

Time delays between the four images of the core of the radio 
triple source, which we believe to be variable, are calculated 
for the mass density profile in Eqn.(l) using the formulae 
of Cooke & Kantowski (1975). The time delay, At, for the 
arrival at the observer of light from an image relative to, say, 
the undeflected source, is the sum of the time delay due to 
differences in geometric path length between the deflected 
and undeflected rays, Ar 9 , and that due to the apparent 
slowing of light by the gravitational potential of the lens, 
At p . With the mass density profile in Eqn.(l), 

^ = ^(f) 2 |I( Zi )| 2 
At, = K36^4HgRe{ f (^X^-l)^ 

In the above expression, the quantity 6< is the value (x 2 + 
y 2 /(l — e 2 ))/a 2 , where x + iy = Zj, the location of im- 
age i relative to the lens centre in the image plane. N = 
(1 + zi)(Di/ Fd), which combination of angular diameter dis- 
tances has dimensions L 1 and calls for additional input: the 
lens and source redshifts and a cosmological model to be 
adopted. The integral in the expression for At p is evaluated 
numerically. Table 3 summarizes the time delay results; for 
the sake of completeness the delays for images arising from 
the lobes of the source are also listed. These are of relevance 
if monitoring reveals relative motion between the core and 
the knots, as is seen in compact symmetric objects (see Sec. 
7.1 ahead). 



6 FITS WITH SUBSETS OF CONSTRAINTS 

It is interesting to ask how the modeling behaves when only 
one set of images is used. As discussed in Section 3.1, the 
positions of the four core images provide by themselves a 



total of 6 constraints, which is exactly equal to the number 
of model parameters. Not surprisingly, an excellent match is 
obtained (4>core = 0.421), with the scale length restricted 
to be less than 10~ 5 asec (effectively a 5 parameter lens). 
The smallness of the mismatch function <f> probably points 
to some degree of redundancy in the information provided 
by images in the core quad; the pairs 1 & 4 and 3 & 6 are 
fairly symmetrically disposed with respect to the lens major 
axis (Table 4). The image magnification ratios 3/1, 4/1 & 
6/1 are 0.95, 1.87 & 0.57 respectively. These should be com- 
pared with the observed values of 0.69, 2.61 & 0.61 for these 
images, recalling that they are known to be highly variable, 
and with the original model values of 1.13, 2.08 & 0.86 re- 
spectively. (The value of 4>core for the core images alone 
in the original model is 3.9 and for all the images together 
is 4>all = 9.8). However, if the scale length is restricted 
and modeling is undertaken of all the images together, the 
fit is a bit worse: 4>all = 15.9 for all the images together, 
or 4>core = 10.3 for the core images alone. Save for the 
scale length in these various models, the other parameters 
are similar; the lens centre shifts by between 10 to 20 mas 
between models, and the model velocity dispersion and ec- 
centricity change by about 4%, which is within the 90% con- 
fidence limits on the original model. The orientation angle 
is almost unaffected. VLB A and MERLIN observations un- 
derway (Sykes et al. 1997) should provide better determined 
constraints from which to build more accurate models in the 
future. 



7 DISCUSSION AND CONCLUSIONS 

Is B1933+503 a lensed system? Despite the absence of opti- 
cal information on the images and a redshift for the source, 
the simplicity with which the lensing picture explains this 
strikingly unusual ten component system indicates that it 
is not just a collection of physically related radio features 
but is indeed a lensed source. This conclusion is underlined 
by the fact that actual lens modeling accounts so well for 
details of the observations. 

The Mass Profile of the Lens Galaxy: The mass model 
adopted (a non-singular isothermal ellipsoid) is only an ap- 
proximation to possibly more complex substructure within 
the lens, especially if it is a disk system as is suggested by 
the high degree of flattening in the optical image. Since three 
source components are lensed by the same galaxy, it will 
provide a strong test of the usual simplifying assumptions 
made when modeling lensed systems. Although a quad sys- 
tem provides a surfeit of constraints (6 positional + 3 image 
flux ratio = 9) over parameters for a simple 6 parameter lens 
model such as the one used in the present work, the images 
typically tend to form at roughly the same radial distance 
from the centre of the lens and convey little information on 
the radial profile or the scale length the mass distribution of 
the galaxy. This is particularly true if the the scale length is 
quite small compared with the scale of the image splitting. 
A double image system manages to sample radially distinct 
portions of the lens, but in general on its own does not pro- 
vide enough constraints (2 positional + 1 image flux ratio 
= 3 quantities) for even the simplest elliptical lens models. 
B 1933+503 is remarkable in having two quads and a double, 
and the fact that the best-fit model calls for a non-singular 
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lens model (scale length of 113 mas, translating to a typical 
value of about 0.5 kpc at zi = 0.755, H = 100 km/s/Mpc) 
is therefore a matter of some significance. 

Predictions for high-resolution observations: The constraints 
used in the present model are independent of the image flux 
density ratios, but modeling relies upon the image parities, 
which are robust constraints. The image parities are usu- 
ally weak constraints in a quad or a double image system, 
but in the present 10-image configuration effectively restrict 
the location of the tangential critical curve. The model in 
this paper makes certain predictions: (i) Magnification ratios 
for the core images are given in Table 2, as obtained from 
modeling. Monitoring variations in the flux density these 
images will provide both the time delays (also predicted by 
the model, see Table 3) and the actual lens magnification 
ratios, (ii) The source of the quad (2a, 2b, 5, 7), a knot 
with a steep radio spectrum, should itself be multicompo- 
nent or extended. In the present model, this source is only 
partially imaged in the vicinity of 2a and 2b. VLB A observa- 
tions at 1.7 GHz have been carried out (Marlow et al. 1997), 
and these should reveal with greater accuracy the correspon- 
dence between the peaks in 2a and 2b and features in the 
other images 5 and 7. 

7.1 The source — a Compact Symmetric Object? 

The source of this remarkable ten-image system appears to 
be a compact radio triple source, an almost linear set of knot, 
core and a second, weaker knot. The stronger and weaker 
knots have separations of 194.5 and 120 mas from the core 
in the present model (see Table 2). This corresponds to a 
physical extent of about 1.3 kpc (H = 100 km/s/Mpc) and 
thus the source could be a compact symmetric object (Read- 
head et al. 1994). Observations over time will indicate, by 
changes in the locations and/or image flux densities of the 
knot images, whether they are moving relative to the cen- 
tral core, as in the case of the radio galaxy 1946+708 (Taylor 
& Vermeulen 1997). While this would no doubt complicate 
the modeling described earlier, it would in itself be a very 
interesting phenomenon to observe. It will be necessary to 
'deconvolve' the effects due to differences in light travel time 
to the observer for the different images from those due to 
the intrinsic source motion. Each epoch of observation will 
catch upto four different phases in the time evolution of the 
source structure in each quad, due to the relative time delays 
between the images. In addition, there is the advantage of 
spatial magnification due to lensing, which can be as high in 
specific directions as factors of 10—20 for images such as 2a, 
2b, and 8. The model in the present paper makes an inter- 
esting prediction in the event that the brighter knot moves 
outward from the core (see Figure 1): the peaks of images 2a 
and 2b should disappear after initially brightening and then 
merging with each other. Image 5 will also move outward 
with respect to image 4. If the fainter knot moves outward 
from the core, image 8 should initially become fainter (it 
moves into a region of lower magnification) ; there will be an 
increase in separation between images 1 and la. 

7.2 B1933+503 as a cosmological probe 

A well-established goal of gravitational lensing is to deter- 
mine values of parameters for cosmological models (Refsdal 



1964). From the expressions in Section the factor H is in- 
versely proportional to Hubble's Constant; thus, with a well- 
fit model for the lens and measured values for the time delays 
between pairs of images, it is possible to estimate Hubble's 
Constant within a specific cosmological model. However, the 
quantity K involves the source redshift, z s , which in the case 
of B1933+503 is unknown. A study of the lens galaxy could 
provide the key to using this system to estimate Hubble's 
Constant. From Section kA it is evident that information in- 
volving the source redshift is restricted to the term Fd- If it 
proves possible to observationally determine the circular ve- 
locity V c of the lens through redshifted HI or molecular line 
studies, then from Equation (4) and the model parameters, 
the quantity Fd — D 3 /Di s is obtainable in principle. This 
is under the assumption that lensing is due to the galaxy 
at zi = 0.755, with no significant contribution from external 
shear or matter along the line of sight. In view of the de- 
gree to which the simple model presented here accounts for 
B1933+503, this does not seem unreasonable. In general, the 
presence of distributed matter along the line of sight cannot 
be ruled out, so at least an upper bound can be put on the 
value of Hubble's Constant. 

The treatment in this paper assumes that the lens is 
viewed edge-on; including an inclination along the line of 
sight is straightforward with regard to the lens model (see 
Bourassa & Kantowski 1975), and should observations pro- 
vide sufficient detail to permit the inclination to be estab- 
lished, B1933+503 will be a rather promising system for 
cosmological studies. 
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Image 


Rel. (RA, Dec) 


Obs. flux ratio at Freq. (GHz) 




(mas) 


1.7 


5 


8.4 


15 


2a 


( 171, 414) 










2b 


( 51, 420) 










5 


(-531, -497) 


1.00 


1.00 


1.00 


1.00 




( ^QK — 1^4*1 
i oyo, — i-O'ij 


1.25 


1.15 


1.16 


1.13 


1 


( 447, 495) 


1.00 


1.00 


1.00 


1.00 


3 


(-389, 158) 


0.69 


0.84 


0.79 


0.61 


4 


(-397, -299) 


2.61 


3.46 


3.65 


3.78 


6 


( 230,-387) 


0.61 


0.96 


0.83 


0.78 


la 


( 545, 584) 


1.00 








8 


(-114, -335) 


4.00 









Table 1. Inputs to the modeling: For each set of images, the positions relative to an initial guess lens centre at (397,299) mas wrt image 
4, and the image flux ratios relative to that image of each set whose value is 1.00, for various frequencies of observation. The model 
uses only the 1.7 GHz data; the other observations are included in the Table to indicate the (un)reliability of various image flux ratios 
as constraints. Values for images 2a & 2b are not listed as it is apparent that these are mappings of only a fraction of the source of 
images 5 & 7 (see text). In the higher frequency observations, the resolution and/or dynamic range was inadequate to pick up 2a & 2b 
as separate components. Images la & 8 are not detected in the higher frequency observations. The core images 1, 3, 4 & 6 show signs of 
variability (c/. discussion in Sykes et al. 1997). Most data are derived from Sykes et al. (1997); positions for 2a & 2b are from Browne 
(private communication) . 



Image Model 


Flux Ratio 


Pm 


Source 


Set Magn. 


(obsvd) 


(-fern) 


Pos. (mas) 


Lobe quad: 








5 2.58 


1.00 


+i 


SRA : 




1.00 


+i 


-105 ± 5 


7 3.05 


1.18 (1.03, 1.62) 


-l 


SDec : 




1.25 (±0.04) 


-l 


-66 ± 1 


2a 




+i 








+i 




2b 




-l 








-l 




Core quad: 








1 2.89 


1.00 


+i 


SRA : 




1.00 


+i 


36 ± 3 


3 3.25 


1.13 (1.02,1.42) 


-l 


SDec : 




0.69 (±0.14) 


-l 


68 ±4 


4 6.01 


2.08 (1.98,2.52) 


+i 






2.61 (±0.31) 


+i 




6 2.48 


0.86 (0.75,1.16) 


-l 






0.61 (±0.13) 


-l 




Lobe double: 








la 2.28 


1.00 


+i 


SRA : 




1.00 


+i 


123 ± 5 


8 14.24 


6.24 (4.62,10.22) 


-l 


SDec : 




4.00 (±1.83) 


-l 


151 ±4 



Table 2. Model Predictions: For the components grouped as discussed in Section 3, the model magnifications, model flux density ratios 
(with 90% confidence intervals determined from a Monte Carlo exercise alongside) — in the succeeding line, these are compared with 
observed values taken from the 1.7 GHz map of Sykes et al. 1997 (with errors as inferred from that paper); the model image parities P m 
(with the theoretically expected values P em on the following line), and the model-recovered source positions (5(RA), <5(Dec)) for each 
image group, w.r.t the lens centre. The source positions are mean values for all backprojected images belonging to a particular group 
(i.e., quad or double). Beside these values are the standard deviations for each set, for the best-fit model itself; ideally, these should be 
as small as possible. 
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Image Time Delay 90% Confidence Interval 

Set (in units K days/Gpc) 

Lobe quad: 

5 0.00 

7 4.32 

2a 3.66 

2b 3.83 



Core quad: 

1 0.00 

3 2.74 

4 2.37 
6 3.15 



Lobe Double: 
la 0.00 
8 6.22 (5.15, 7.02) 



Table 3. Predicted time delays for the three image sets: The image listed at top in each set is the reference image corresponding to the 
global minimum in arrival times for light from the source (core or a lobe). H = (1 + zi)(Di/F^), as discussed in Section 5 (see text); a 
typical value is 3 Gpc). The last column lists the 90% confidence intervals for each value, from the Monte Carlo exercise. 



Image B value 


\4>\ 


Set (asec) 


(degrees) 


Lobe Quad: 




5 1.24 


83 


7 0.46 


27 


2a 0.75 


65 


2b 0.66 


51 


Core Quad: 




1 1.14 


86 


3 0.50 


19 


4 0.80 


75 


6 0.44 


17 


Lobe Double: 




la 0.80 


3 


8 0.37 


19 



Table 4. Illustrating the extent of sampling of the lens plane by the various image systems. B is the elliptical radius of the given image: 
in Cartesian coordinates, referred to the lens centre and with the x— axis along the lens major axis, B 2 = x 2 + j/ 2 /(l — e 2 ). \ ip , the 
norm of the angle each image makes with respect to the lens major axis (when referred to the lens centre), is listed alongside. Images of 
similar | tp | and B values provide similar constraints on the modeling if they belong to the same image system. 



(3.56, 4.89) 
(3.14, 4.05) 
(3.20, 4.29) 



(2.11, 3.13) 
(1.95, 2.71) 
(2.52, 3.58) 
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Figure 1. Lens model for B1933+503. (a) : Image plane critical curves for the lens model, on which observed image positions have been 
superposed. The lens is centred at (0,0). Images landing very near a critical curve suffer very high magnifications, tending to merge with 
nearby images across such curves (2a with 2b, for example). Image positions marked by open circles with central dots form the quad 
arising from lensing of the source's core, open circles mark one lobe's quad images, and the small filled squares denote the other lobe's 
double images, (b) : The corresponding source plane caustics for the model. The positions of the three source components are marked 
with symbols corresponding to their images, a circle with a dot for the core of the source, an open circle for the stronger lobe and a small 
filled square for the weaker one. The cross superposed on the open circle is the source location when the positions of images 2a and 2b 
are included (see Section 3.1); the slight offset with respect to the original source position, obtained from images 5 and 7 alone, arises 
because images 2a and 2b represent only a fraction of the source that gets mapped onto 5 and 7. 
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Figure 2. Results from the Monte Carlo exercise (Section 4). The plots are histograms of binned values of the six lens parameters, and 
of the mismatch function rj>. The plot for the mismatch function shows the buildup of the distribution as the number of experiments, 
N eX p, is increased from 10 3 to 10 4 . Note that the lens parameters are varied simultaneously in this set of experiments, which results in 
a significantly wider spread in recovered values than is currently recognised in the literature. A spread of the order of 30 % is apparent 
in the predicted time delay values as well (Table 3). 
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